2.1 One Dimension Summary of Data

• 5-number summary: summary(data) = returns min, 1st quartile, median, mean, 3rd quartile, max
• Boxplots: boxplot(data, col = “blue”) = produces a box with middles 50% highlighted in the specified color
  • whiskers = \(\pm 1.58IQR/\sqrt{n}\)
    • IQR = interquartile range, Q\(_3\) - Q\(_1\)
  • box = 25%, median, 75%
• Histograms: histograms(data, col = “green”) = produces a histogram with specified breaks and color
  • breaks = 100 = the higher the number is the smaller/narrower the histogram columns are
• Density plot: rug(data) = density plot, add a strip under the histogram indicating location of each data point

• Barplot: barplot(data, col = wheat) = produces a bar graph, usually for categorical data

• Overlaying Features
  • abline(h/v = 12) = overlays horizontal/vertical line at specified location
    • col = “red” = specifies color
    • lwd = 4 = line width
    • lty = 2 = line type
  • abline(v = median(data), col = "magenta", lwd = 4)

2.2 Two Dimensional Summaries

• multiple/overlay 1D plots (using lattice/ggplot2)
• Multiple boxplots:

boxplot(pm25 ~ region, data = pollution, col = “red”)

• Multiple Histogram:

par(mfrow = c(2, 1), mar = c(4, 4, 2, 1))                       # = set margin
hist(subset(pollution, region == "east")$pm25, col = "green")   # = first histogram
hist(subset(pollution, region == "west")$pm25, col = "green")   # = second histogram

• Scatterplot

with(pollution, plot(latitude, pm25, col = region))
abline(h = 12, lwd = 2, lty = 2)        # = plots horizontal dotted line
plot(jitter(child, 4)~parent, galton)   # = spreads out data points at the same position to simulate 
                                        #   measurement error/make high frequency more visibble

• Multiple scatterplots

par(mfrow = c(1, 2), mar = c(5, 4, 2, 1))                                           # = sets margins
with(subset(pollution, region == "west"), plot(latitude, pm25, main = "West"))      # = left scatterplot
with(subset(pollution, region == "east"), plot(latitude,  pm25, main = "East"))     # = right scatterplot

2.3 Process of Making a Plot/Considerations

• where will plot be made? screen or file?
• how will plot be used? viewing on screen/web browser/print/presentation?
• large amount of data vs few points?
• need to be able to dynamically resize?
• What graphics system to use? 3 core plotting systems in R:
  • base: “artist’s palette” model
  • lattice: Entire plot specified by one function; conditioning
  • ggplot2: Mixes elements of Base and Lattice

3.1 Base Graphics Functions and Parameters

• Arguments
  • lty: line type (default is solid)
  • pch: plotting symbol (default = open circle)
    • 0=blank, 1=solid (default), 2=dashed, 3=dotted, 4=dotdash, 5=longdash, 6=twodash
  • lwd: line width (integer)
  • col: plotting color (number string or hexcode, colors() returns vector of colors)
  • xlab, ylab: x-y label character strings
  • cex: numerical value giving the amount by which plotting text/symbols should be magnified relative to the default
    • cex = 0.15 * variable: plot size as an additional variable
• par() function = specifies global graphics parameters, affects all plots in an R session (can be overridden)
  • las: orientation of axis labels
  • bg: background color
  • mar: margin size (order = bottom left top right)
  • oma: outer margin size (default = 0 for all sides)
  • mfrow: number of plots per row, column (plots are filled row-wise)
  • mfcol: number of plots per row, column (plots are filled column-wise)
  • can verify all above parameters by calling par("parameter")
  • list of all parameters by calling names(par())
  • use dev.off or plot.new to reset to the defaults
• plotting functions
  • lines: adds lines to a plot, given a vector of x values and corresponding vector of y values
  • points: adds a point to the plot
  • text: add text labels to a plot using specified x,y coordinates
  • title: add annotations to x,y axis labels, title, subtitles, outer margin
  • mtext: add arbitrary text to margins (inner or outer) of plot
  • axis: specify axis ticks
  • legend: explains what the symbols mean
  • abline: add a horizontal/vertical line

3.2 Base Plot Examples

• Basic Examples:

Scatterplot

library(datasets)
data(cars)
with(cars, plot(speed, dist))

Histogram

library(datasets)
hist(airquality$Ozone)

Boxplot

library(datasets)
airquality <- transform(airquality, Month = factor(Month))
boxplot(Ozone ~ Month, airquality, xlab = "Month", ylab = "Ozone (ppb)")

• Advanced Examples:

Scatterplot with Regression Line

library(datasets)
# type =“n” sets up the plot and does not fill it with data
with(airquality, plot(Wind, Ozone, main = "Ozone and Wind in New York City", type = "n"))

# subsets of data are plotted here using different colors
with(subset(airquality, Month == 5), points(Wind, Ozone, col = "blue"))
with(subset(airquality, Month != 5), points(Wind, Ozone, col = "red"))
legend("topright", pch = 1, col = c("blue", "red"), legend = c("May", "Other Months"))

# regression line is produced here
model <- lm(Ozone ~ Wind, airquality)
abline(model, lwd = 2)

Scatterplot with factor variables

x <- rnorm(100)
y <- x + rnorm(100)
g <- gl(2, 50, labels = c("Male", "Female"))

plot(x, y, type = "n")
points(x[g == "Male"], y[g == "Male"], col = "green")
points(x[g == "Female"], y[g == "Female"], col = "blue", pch = 19)

3.3 Multiple Plot Example

• Note: typing example(points) in R will launch a demo of base plotting system and may provide some helpful tips on graphing

# this expression sets up a plot with 1 row 3 columns, sets the margin and outer margins
par(mfrow = c(1, 3), mar = c(4, 4, 2, 1), oma = c(0, 0, 2, 0))
with(airquality, {
    # here three plots are filled in with their respective titles
    plot(Wind, Ozone, main = "Ozone and Wind")
    plot(Solar.R, Ozone, main = "Ozone and Solar Radiation")
    plot(Temp, Ozone, main = "Ozone and Temperature")
    # this adds a line of text in the outer margin*
    mtext("Ozone and Weather in New York City", outer = TRUE)}
)

5.1 lattice Functions and Parameters

• Funtions
  • xyplot() = main function for creating scatterplots
  • bwplot() = box and whiskers plots (box plots)
  • histogram() = histograms
  • stripplot() = box plot with actual points
  • dotplot() = plot dots on “violin strings”
  • splom() = scatterplot matrix (like pairs() in base plotting system)
  • levelplot()/contourplot() = plotting image data
• Arguments for xyplot(y ~ x | f * g, data, layout, panel)
  • default blue open circles for data points
  • formula notation is used here (~) = LHS is the y-axis variable, and the RHS is the x-axis variable
  • f/g = conditioning/categorical variables (optional)
    • basically creates multi-panelled plots (for different factor levels)
    • * indicates interaction between two variables
    • intuitively, the xyplot displays a graph between x and y for every level of f and g
    • Example: xyplot(price~carat | color*cut, data = diamonds, strip = FALSE, pch = 20)
  • data = the data frame/list from which the variables should be looked up
    • if nothing is passed, the parent frame is used (searching for variables in the workspace)
    • if no other arguments are passed, defaults will be used
  • layout = specifies how the different plots will appear
    • layout = c(5, 1) = produces 5 subplots in a horizontal fashion
    • padding/spacing/margin automatically set
  • [optional] panel function can be added to control what is plotted inside each panel of the plot
    • panel functions receive x/y coordinates of the data points in their panel (along with any additional arguments)
    • ?panel.xyplot = brings up documentation for the panel functions
    • panel.abline
    • panel.lmline
  • some arguments are same as in base plotting: xyplot(Ozone ~ Wind, data = airquality, pch=8, col="red", main="Big Apple Data")
  • Note: no base plot functions can be used for lattice plots

5.2 lattice Examples

• Basic Example:

library(lattice)
state <- data.frame(state.x77, region = state.region)
xyplot(Life.Exp ~ Income | region, data = state, layout = c(4,1))

• Advanced Example (Lattice Panel Functions):

library(lattice)
set.seed(10)
x <- rnorm(100)
f <- rep(0:1, each = 50)
y <- x + f - f * x+ rnorm(100, sd = 0.5)
f <- factor(f, labels = c("Group 1", "Group 2"))
# Plot with 2 panels with custom panel function
xyplot(y ~ x | f, panel = function(x, y, ...) {
    # call the default panel function for xyplot
    panel.xyplot(x, y, ...)
    # adds a horizontal line at the median
    panel.abline(h = median(y), lty = 2)
    # overlays a simple linear regression line
    panel.lmline(x, y, col = 2)
})

6.1 ggplot2 Functions and Parameters

• basic components of a ggplot2 graphic
  • data frame = source of data
  • aesthetic mappings = how data are mappped to color/size (x vs y)
  • geoms = geometric objects like points/lines/shapes to put on page
  • facets = conditional plots using factor variables/multiple panels
  • stats = statistical transformations like binning/quantiles/smoothing
  • scales = scale aesthetic map uses (i.e. male = red, female = blue)
  • coordinate system = system in which data are plotted

library(ggplot2)
qplot(displ, hwy, data = mpg, color = drv, shape = drv)

qplot(y = hwy, data = mpg, color = drv)

qplot(displ, hwy, data = mpg, geom = c("point", "smooth"), method="loess")

load("data/maacs.Rda")
qplot(log(pm25), log(eno), data = maacs, color = mopos) + geom_smooth(method = "lm")

qplot(drv,hwy, data = mpg, geom = "boxplot")

# Can even be subdivided into several boxes:
qplot(drv,hwy, data = mpg, geom = "boxplot", color = manufacturer)

qplot(hwy, data = mpg, fill = drv)

library(gridExtra)
q1 <- qplot(log(eno), data = maacs, geom="density")
q2 <- qplot(log(eno), data = maacs, geom="density", color = mopos)
grid.arrange(q1, q2, ncol = 2)

qplot(displ, hwy, data = mpg, facets = . ~ drv)

qplot(hwy, data = mpg, facets = drv ~ ., binwidth = 2)

library(datasets)
data(airquality)
# Need to transform 'Month' to a factor variable first
airquality = transform(airquality, Month = factor(Month))
qplot(Wind, Ozone, data = airquality, facets = . ~ Month)

6.2 ggplot2 Examples

• Basic Example (qplot):

library(ggplot2)
data(mpg)
qplot(displ, hwy, data = mpg)

• Comprehensive Example (ggplot):

# initiates ggplot
g <- ggplot(maacs, aes(logpm25, NocturnalSympt))
g + geom_point(alpha = 1/3)                             # adds points (alpha -> transparency)
      - facet_wrap(bmicat ~ no2dec, nrow = 2, ncol = 4)     # make panels
      - geom_smooth(method="lm", se=FALSE, col="steelblue") # adds smoother
      - theme_bw(base_family = "Avenir", base_size = 10)    # change theme
      - labs(x = expression("log " * PM[2.5])               # add labels
      - labs(y = "Nocturnal Symptoms”)
      - labs(title = "MAACS Cohort”)

7.1 Procedure for Constructing Hierarchical Clusters (hclust function)

1. calculate all pair wise distances between all points to see which points are closest together
   • dist(data.frame(x=x, y=y) = returns pair wise distances for all of the (x,y) coordinates
   • Note: dist() function uses Euclidean distance by default
2. group two closest points from the calculated distances and merge them to a single point
3. find the next two closest points and merge them, and repeat
4. order of clustering is shown in the dendrogram

7.2 Approaches for Merging Points/Clusters

• the approach is specified in the argument method = "complete" or "average" in hclust() function
• average linkage = taking average of the x and y coordinates for both points/clusters (center of mass effectively)
• complete linkage = to measure distance of two clusters, take the two points in the clusters that are the furthest apart
• Note: two approaches may produce different results so it’s a good idea to use both approaches to validate results

7.3 Characteristics of Hierarchical Clustering Algorithms

• clustering result/plot maybe unstable
  • changing few points/outliers could lead to large changes
  • change different distance metrics to see how sensitive the clustering is
  • change merging strategy
  • scaling of variables could affect the clustering (if one unit/measurement is much larger than another)
• deterministic = running the hclust function with same parameters and the same data will produce the same plot
• determining how many clusters there are (where to cut) may not always be clear
• primarily used for exploratory data analysis, to see over all pattern in data if there is any at all

7.4 hclust Function and Example

• hh <- hclust(dist(dataFrame)) function = produces a hierarchical cluster object based on pair wise distances from a data frame of x and y values
  • dist() = defaults to Euclidean, calculates the distance/similarity between two observations; when applied to a data frame, the function applies the \(\sqrt{(A_1 - A_2)^2 + (B_1 - B_2)^2 + ... + (Z_1 -Z_2)^2 }\) formula to every pair of rows of data to construct a matrix of distances between the rows
    • order of the hierarchical cluster is derived from the distance
  • plot(hh) = plots the dendrogram
  • automatically sorts column and row according to cluster
  • names(hh) = returns all parameters of the hclust object
    • hh$order = returns the order of the rows/clusters from the dendrogram
    • hh$dist.method = returns method for calculating distance/similarity
• Note: dendrogram that gets generated DOES NOT show how many clusters there are, so cutting (at 2.0 level for example) must be done to determine number of clusters — must be a convenient and sensible point

# create points and plot them
set.seed(1234)
par(mar = c(1, 1, 1, 2), mfrow=c(1,2))
x <- rnorm(12, mean = rep(1:3, each = 4), sd = 0.2)
y <- rnorm(12, mean = rep(c(1, 2, 1), each = 4), sd = 0.2)
plot(x, y, col = "blue", pch = 19, cex = 2)
text(x + 0.05, y + 0.05, labels = as.character(1:12))

# construct hierarchical clusters and plot dendrogram
dataFrame <- data.frame(x = x, y = y)
distxy <- dist(dataFrame)               # calculate all pairwise distances
hClustering <- hclust(distxy)
plot(hClustering)

7.5 myplcclust Function and Example

• Note: myplcclust = a function to plot hclust objects in color (clusters labeled 1 2 3 etc.), but must know how many clusters there are initially

myplclust <- function(hclust, lab = hclust$labels,
    lab.col = rep(1, length(hclust$labels)), hang = 0.1, ...) {
    # modifiction of plclust for plotting hclust objects *in colour*! Copyright
    # Eva KF Chan 2009 Arguments: hclust: hclust object lab: a character vector
    # of labels of the leaves of the tree lab.col: colour for the labels;
    # NA=default device foreground colour hang: as in hclust & plclust Side
    # effect: A display of hierarchical cluster with coloured leaf labels.
    y <- rep(hclust$height, 2)
    x <- as.numeric(hclust$merge)
    y <- y[which(x < 0)]
    x <- x[which(x < 0)]
    x <- abs(x)
    y <- y[order(x)]
    x <- x[order(x)]
    plot(hclust, labels = FALSE, hang = hang, ...)
    text(x = x, y = y[hclust$order] - (max(hclust$height) * hang), labels = lab[hclust$order],
 col = lab.col[hclust$order], srt = 90, adj = c(1, 0.5), xpd = NA, ...)
}
# example
dataFrame <- data.frame(x = x, y = y)
distxy <- dist(dataFrame)
hClustering <- hclust(distxy)
myplclust(hClustering, lab = rep(1:3, each = 4), lab.col = rep(1:3, each = 4))

7.6 heatmap Function and Example

• tutorial on heatmaps
• heatmap(data.matrix) function = similar to image(t(x))
  • good for visualizing high-dimension matrix data, runs hierarchical analysis on rows and columns of table
  • yellow = high value, red = low value
  • Note: the input must be a numeric matrix, so as.matrix(data.frame) can be used to convert if necessary

dataFrame <- data.frame(x = x, y = y)
set.seed(143)
dataMatrix <- as.matrix(dataFrame)[sample(1:12), ]
# heatmap(dataMatrix, col=cm.colors(19))
heatmap(dataMatrix)

7.7 image Function and Example

• image(x, y, t(dataMatrix)[, nrow(dataMatrix):1]) = produces similar color grid plot as the heatmap() without the dendrograms
  • t(dataMatrix)[, nrow(dataMatrix)]
    • t(dataMatrix) = transpose of dataMatrix, this is such that the plot will be displayed in the same fashion as the matrix (rows as values on the y axis and columns as values on the x axis)
      • example 40 x 10 matrix will graph the 10 columns as x values and 40 rows as y values
    • [, nrow(dataMatrix)] = subsets the data frame in reverse column order; when combined with the t() function, it reorders the rows of data from 40 to 1, such that the data from the matrix is displayed in order from top to bottom
      • Note: without this statement the rows will be displayed in order from bottom to top, as that is in line with the positive y axis
  • x, y = used to specify the values displayed on the x and y axis
    • Note: must be in increasing order

image(1:2, 1:12, t(dataMatrix)[, nrow(dataMatrix):1])

8.1 Procedure for Constructing K-means Clusters (kmeans function)

1. choose three random points as the starting centroids
2. take each of the data points and assign it to the closest centroid (creating a cluster around each starting point)
3. take each cluster and recalculate the centroid (taking the mean) with its enclosed data points
4. repeat step 2 and 3 (reassign points to centroids and update centroid locations) until a stable result is achieved

set.seed(1234)
x <- rnorm(12,mean=rep(1:3,each=4),sd=0.2)
y <- rnorm(12,mean=rep(c(1,2,1),each=4),sd=0.2)
dataFrame <- data.frame(x=x,y=y)

# specifies initial number of clusters to be 3
kmeansObj <- kmeans(dataFrame,centers=3)
names(kmeansObj)

## [1] "cluster"      "centers"      "totss"        "withinss"    
## [5] "tot.withinss" "betweenss"    "size"         "iter"        
## [9] "ifault"

# returns cluster assignments
kmeansObj$cluster

##  [1] 3 3 3 3 1 1 1 1 2 2 2 2

# plotting
par(mar=rep(0.2,4))
plot(x,y,col=kmeansObj$cluster,pch=19,cex=2)
points(kmeansObj$centers,col=1:3,pch=3,cex=3,lwd=3)

8.2 K-means clustering and heatmaps

• similar to hierarchical clustering

set.seed(1234)
dataMatrix <- as.matrix(dataFrame)[sample(1:12),]
kmeansObj <- kmeans(dataMatrix,centers=3)

# plotting
par(mfrow=c(1,2), mar = c(2, 4, 0.1, 0.1))
image(t(dataMatrix)[,nrow(dataMatrix):1],yaxt="n")
image(t(dataMatrix)[,order(kmeansObj$cluster)],yaxt="n")

8.3 Characteristics of K-means Clustering Algorithms

• requires number of clusters initially
  • pick by eye/intuition
  • pick by cross validation/information theory, etc.
• not deterministic (starting points chosen at random)
  • useful to run the algorithms a few times with different starting points to validate results

9.1 Singular Value Decomposition (SVD)

• Let \(X\) = matrix which each variable in column (measurement) and each observation in row (subject)
• SVD in this case is a matrix decomposition process, in which X is divided into three separate matrices as follows: \[X = UDV^T\]
  • \(U\) = left singular vector, orthogonal matrix (columns independent of each other)
  • \(D\) = singular values, diagonal matrix
  • \(V\) = right singular vector, orthogonal matrix (columns independent of each other)
  • Note: orthogonal implies that a matrix is always invertible [\(A^{-1} = A^T\)] and that the product of the matrix and its transpose equals the identity matrix [\(AA^T = I\)]
    • when a orthogonal matrices, \(A\), is multiplied by another matrix, \(B\), it is effectively a linear transformation in that the length and angles of \(B\) are preserved
  • Note: diagonal implies that any value outside of the main diagonal (\(\searrow\)) = 0 \[A = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 2 & 0 \\ 0 & 0 & 3 \end{bmatrix}\]
• Note: scale of data matters for SVD/PCA (scaling the data may help), patterns detected maybe mixed together, and computation is intensive for these operations

9.2 Principal Components Analysis (PCA)

• first scale the variables and run SVD on normalized matrix
  • scaling = subtract each column by its mean and divide by its standard deviation
• principal components = the right singular values or the \(V\) matrix

9.3 SVD and PCA Example

• \(U\) and \(V\) Matrices
  • s <- svd(data) = performs SVD on data (\(n \times m\) matrix) and splits it into \(u\), \(v\), and \(d\) matrices
    • s$u = \(n \times m\) matrix \(\rightarrow\) horizontal variation
    • s$d = \(1 \times m\) vector \(\rightarrow\) vector of the singular/diagonal values
      • diag(s$d) = \(m \times m\) diagonal matrix
    • s$v = \(m \times m\) matrix \(\rightarrow\) vertical variation
    • s$u %*% diag(s$d) %*% t(s$v) = returns the original data \(\rightarrow\) \(X = UDV^T\)
  • scale(data) = scales the original data by subtracting each data point by its column mean and dividing by its column standard deviation

# running svd
svd1 <- svd(scale(dataOrdered))

# create 1 by 3 panel plot
par(mfrow=c(1,3))
# data heatmap (sorted)
image(t(dataOrdered)[,nrow(dataOrdered):1])
# U Matrix - first column
plot(svd1$u[,1],40:1,,xlab="Row",ylab="First left singular vector",pch=19)
# V vector - first column
plot(svd1$v[,1],xlab="Column",ylab="First right singular vector",pch=19)

• \(D\) Matrix and Variance Explained
  • \(d\) matrix (s$d vector) captures the singular values, or variation in data that is explained by that particular component (variable/column/dimension)
  • Recall that \(D\) is the diagonal matrix sandwiched in between \(U\) and \(V^T\) in the SVD representation of the data matrix. The diagonal entries of \(D\) are like weights for the \(U\) and \(V\) columns accounting for the variation in the data. They’re given in decreasing order from highest to lowest.
  • proportion of variance Explained = converting the singular values to variance (square the values) and divide by the total variance (sum of the squared singular values)
    • effectively the same pattern as the singular values, just converted to percentage
    • in this case, the first component/dimension, which captures the shift in means (see previous plot) of SVD captures about 40% of the variation

# create 1 x 2 panel plot
par(mfrow=c(1,2))
# plot singular values
plot(svd1$d,xlab="Column",ylab="Singular value",pch=19)
# plot proportion of variance explained
plot(svd1$d^2/sum(svd1$d^2),xlab="Column",ylab="Prop. of variance explained",pch=19)

Variance explained by single variable:

constantMatrix <- dataOrdered*0
for(i in 1:dim(dataOrdered)[1]){constantMatrix[i,] <- rep(c(0,1),each=5)}
svd1 <- svd(constantMatrix)

par(mfrow=c(1,3))
image(t(constantMatrix)[,nrow(constantMatrix):1])
plot(svd1$d,xlab="Column",ylab="Singular value",pch=19)
plot(svd1$d^2/sum(svd1$d^2),xlab="Column",ylab="Prop. of variance explained",pch=19)

• Relationship to PCA
  • p <- prcomp(data, scale = TRUE) = performs PCA on data specified
    • scale = TRUE = scales the data before performing PCA
    • returns prcomp object
    • summary(p) = prints out the principal component’s standard deviation, proportion of variance, and cumulative proportion
  • PCA’s rotation vectors are equivalent to their counterparts in the V matrix from the SVD

# SVD
svd1 <- svd(scale(dataOrdered))

# PCA
pca1 <- prcomp(dataOrdered,scale=TRUE)
# alternatively: pca1 <- prcomp(scale(dataOrdered))

# Plot the rotation from PCA (Principal Components) vs v vector from SVD
plot(pca1$rotation[,1],svd1$v[,1],pch=19,xlab="Principal Component 1",
    ylab="Right Singular Vector 1")
abline(c(0,1))

# summarize PCA
summary(pca1)

## Importance of components%s:
##                          PC1    PC2    PC3    PC4     PC5     PC6     PC7
## Standard deviation     1.995 1.2458 1.0781 1.0091 0.93835 0.72087 0.62791
## Proportion of Variance 0.398 0.1552 0.1162 0.1018 0.08805 0.05197 0.03943
## Cumulative Proportion  0.398 0.5532 0.6694 0.7712 0.85926 0.91122 0.95065
##                            PC8     PC9    PC10
## Standard deviation     0.47621 0.38454 0.34473
## Proportion of Variance 0.02268 0.01479 0.01188
## Cumulative Proportion  0.97333 0.98812 1.00000

• More Complex Patterns
  • SVD can be used to detect unknown patterns within the data (we rarely know the true distribution/pattern about the population we’re analyzing)
  • however, it may be hard to pinpoint exact patterns as the principal components may confound each other
    • in the example below, you can see that the two principal components that capture the most variation have both horizontal shifts and alternating patterns captured in them

Introducing more complex patterns:

set.seed(678910)
# setting pattern
data <- matrix(rnorm(400), nrow = 40)
for(i in 1:40){
  # flip a coin
  coinFlip1 <- rbinom(1,size=1,prob=0.5)
  coinFlip2 <- rbinom(1,size=1,prob=0.5)

  # if coin is heads add a common pattern to that row
  if(coinFlip1){
    data[i,] <- data[i,] + rep(c(0,5),each=5)
  }
  if(coinFlip2){
    data[i,] <- data[i,] + rep(c(0,5),5)
  }
}
hh <- hclust(dist(data))
dataOrdered <- data[hh$order, ]

# perform SVD
svd2 <- svd(scale(dataOrdered))

# true patterns
par(mfrow=c(2,3))
image(t(dataOrdered)[,nrow(dataOrdered):1])
plot(rep(c(0,1),each=5),pch=19,xlab="Column", main="True Pattern 1")
plot(rep(c(0,1),5),pch=19,xlab="Column",main="True Pattern 2")

# v and patterns of variance in rows
image(t(dataOrdered)[,nrow(dataOrdered):1])
plot(svd2$v[,1],pch=19,xlab="Column",ylab="First right singular vector", main="Detected Pattern 1")
plot(svd2$v[,2],pch=19,xlab="Column",ylab="Second right singular vector", main="Detected Pattern 2")

# d and variance explained
svd1 <- svd(scale(dataOrdered))
par(mfrow = c(1, 2))
plot(svd1$d, xlab = "Column", ylab = "Singular value", pch = 19)
plot(svd1$d^2/sum(svd1$d^2), xlab = "Column", ylab = "Percent of variance explained", pch = 19)

• Missing Data
  • SVD cannot be performed on dataset with NA values
  • impute package from Bioconductor can help approximate missing values from surrounding values
    • impute.knn function takes the missing row and imputes the data using the k nearest neighbors to that row
      • k=10 = default value (take the nearest 10 rows)

library(impute)  ## Available from http://bioconductor.org
data2 <- dataOrdered

# set random samples = NA
data2[sample(1:100,size=40,replace=FALSE)] <- NA
data2 <- impute.knn(data2)$data
svd1 <- svd(scale(dataOrdered))
svd2 <- svd(scale(data2))

# plotting
par(mfrow=c(1,2))
plot(svd1$v[,1],pch=19, main="Original")
plot(svd2$v[,1],pch=19, main="Imputed")

9.4 Create Approximations/Data Compression

• SVD can be used to create lower rank representation, or compressed representation of data
• if we look at the variance explained plot below, most of the variation is explained by the first few principal components

# load faceData
load("data/face.rda")

# perform SVD
svd3 <- svd(scale(faceData))
plot(svd3$d^2/sum(svd3$d^2),pch=19,xlab="Singular vector",ylab="Variance explained")

• approximations can thus be created by taking the first few components and using matrix multiplication with the corresponding \(U\), \(V\), and \(D\) components

# Here svd3$d[1] is a constant
approx1 <- svd3$u[,1] %*% t(svd3$v[,1]) * svd3$d[1]
# In these examples we need to make the diagonal matrix out of d
approx5 <- svd3$u[,1:5] %*% diag(svd3$d[1:5])%*% t(svd3$v[,1:5])
approx10 <- svd3$u[,1:10] %*% diag(svd3$d[1:10])%*% t(svd3$v[,1:10])

# create 1 x 4 panel plot
par(mfrow=c(1,4))

# plot original facedata
image(t(approx1)[,nrow(approx1):1], main = "1 Component")
image(t(approx5)[,nrow(approx5):1], main = "5 Component")
image(t(approx10)[,nrow(approx10):1], main = "10 Component")
image(t(faceData)[,nrow(faceData):1], main = "Original")

9.5 Notes and further resources

• Scale matters: when reducing dimensions you have to pay attention to the scales on which different variables are measured and make sure that all your data is in consistent units
• PC’s/SV’s may mix real patterns (see example of more complex patterns, True vs detected patterns)
• Can be computationally intensive
• Tutorial on PCS
• Advanced data analysis from an elementary point of view
• Elements of statistical learning
• Alternatives:
  • Factor analysis
  • Independent components analysis
  • Latent semantic analysis

10.1 grDevices Package

• colors() function = lists names of colors available in any plotting function
• grDevices package has two functions that take palettes of colors and help to interpolate between the colors
  • colorRamp()
  • colorRampPalette()
• The colorRamp and colorRampPalette functions can be used in conjunction with color palettes (eg. RColorBrewer package) to connect data to colors
• colorRamp function
  • takes any set of colors and return a function that takes values between 0 and 1, indicating the extremes of the color palette (e.g. see the gray function)
  • pal <- colorRamp(c("red", "blue")) = defines a colorRamp function
  • pal(0) returns a 1x3 matrix containing values for RED, GREEN, and BLUE values that range from 0 to 255
  • pal(seq(0, 1, len = 10)) returns a 10x3 matrix of 10 colors that range from RED to BLUE (two ends of spectrum defined in the object)

# define colorRamp function
pal <- colorRamp(c("red", "blue"))
# create a color
pal(0.67)

##       [,1] [,2]   [,3]
## [1,] 84.15    0 170.85

• colorRampPalette function
  • takes any set of colors and return a function that takes integer arguments and returns a vector of colors interpolating the palette (like heat.colors or topo.colors)
  • pal <- colorRampPalette(c("red", "yellow")) defines a colorRampPalette function
  • pal(10) returns 10 interpolated colors in hexadecimal format that range between the defined ends of spectrum

# define colorRampPalette function
pal <- colorRampPalette(c("red", "yellow"))
# create 10 colors
pal(10)

##  [1] "#FF0000" "#FF1C00" "#FF3800" "#FF5500" "#FF7100" "#FF8D00" "#FFAA00"
##  [8] "#FFC600" "#FFE200" "#FFFF00"

• rgb function
  • The rgb function can be used to produce any color via red, green, blue proportions, while color transparency can be added via the alpha parameter to rgb
    • red, green, and blue arguments = values between 0 and 1
    • alpha = 0.5 = transparency control, values between 0 and 1
  • returns hexadecimal string for color that can be used in plot/image commands
  • colorspace package can be used for different control over colors

x <- rnorm(500); y <- rnorm(500)
par(mfrow=c(1,2))
# normal scatter plot
plot(x, y, pch = 19, main = "Default")
# using transparency shows data much better
plot(x, y, col = rgb(0, 0, 0, 0.2), pch = 19, main = "With Transparency")

10.2 RColorBrewer Package

• can be found on CRAN that has predefined color palettes
  • library(RColorBrewer)
• 3 types of palettes
  • Sequential = numerical/continuous data that is ordered from low to high
  • Qualitative = categorical data/factor variables
  • Diverging = data that deviate from a value, increasing in two directions (i.e. standard deviations from the mean)
• palette information from the RColorBrewer package can be used by colorRamp() and colorRampPalette() functions
• available colors palettes:

• brewer.pal(n, "BuGn") function
  • n = number of colors to generated
  • "BuGn" = name of palette
    • ?brewer.pal list all available palettes to use
  • returns list of n hexadecimal colors

library(RColorBrewer)
# generate 3 colors using brewer.pal function
cols <- brewer.pal(3, "BuGn")       # returns a vector of colors in hexadecimal format
pal <- colorRampPalette(cols)

par(mfrow=c(1,3))
# heat.colors/default
image(volcano, main = "Heat.colors/Default")
# topographical colors
image(volcano, col = topo.colors(20), main = "Topographical Colors")
# RColorBrewer colors
image(volcano, col = pal(20), main = "RColorBrewer Colors")

• smoothScatter function
  • used to plot large quantities of data points
  • creates 2D histogram of points and plots the histogram
  • default color scheme = “Blues” palette from RColorBrewer package

x <- rnorm(10000); y <- rnorm(10000)
smoothScatter(x, y)